#include "vainshtein.h"
#include <stdio.h>
#include <gsl/gsl_linalg.h>

#define NEAR(a,ipart) (modf(a,&ipart)>0.5)?ceil(a):floor(a)
double ipart;

//void zerosCoax(double *,int,int);
void viewinfo (QTextStream *, int, gsl_matrix_complex *, gsl_vector_complex *, gsl_vector_complex *, gsl_permutation *);

void Vainshtein::stepByStep(QTextBrowser *stack)
{//! This function for approx
    int N,sign=1,itM,itH;   // iterators
    QString str,str1;
    QFile writefile("ANALITIC.CSV");
    QTextStream stream(&writefile);
    //initialisation
    difQ[2]=difQ[1]*(difQ[3]-1);//! convergence
    N=difQ[1]+difQ[2];    //without TEM
    for (int j=80; j<=difQ[2]; j++) HMN[j-1]=M_PI*j/(difQ[3]-1); //HMN[j-1]=zerosCoaxH(mI,j);
    for (int j=80; j<=N; j++) MMN[j-1]=M_PI*(j+0.25); //MMN[j-1]=zerosJs(mI,j);
    gsl_complex z,Aq,beta1,beta,delta;
    gsl_vector_complex * Xnew = gsl_vector_complex_alloc (N);
    //eps=gsl_complex_mul_real(eps,-pow(k,2));
    beta1=gsl_complex_sqrt_real(-pow(k,2)+pow(MMN[nI-1]/a,2));
    beta1=gsl_complex_mul_real(beta1,-1);
    if(writefile.open(QIODevice::WriteOnly))
    {
        for(int p=0; p<difQ[1]; p++)
        {   // for each wave in first area
            itM=0; itH=0;
            beta=gsl_complex_sqrt_real(-pow(k,2)+pow(MMN[p]/a,2));
            //GSL_SET_COMPLEX(&Aq,gsl_sf_bessel_Jn(2,MMN[0])/gsl_sf_bessel_Jn(2,MMN[p]),0);
            if(GSL_REAL(beta)==0) beta=gsl_complex_mul_real(beta,-1);
            GSL_SET_COMPLEX(&z,((a*log(a/(Shield-a))-Shield*log(Shield/(Shield-a)))-a/(p+1))/M_PI,0);
            z=gsl_complex_mul(z,gsl_complex_add(beta1,beta));
            Aq=gsl_complex_mul_real(gsl_complex_exp(z),gsl_sf_bessel_Jn(2,MMN[nI-1])/gsl_sf_bessel_Jn(2,MMN[p]));
            stream << "p \t exp \t div \t Aq \n";
            for(int j=0; j<difQ[1]; j++)
            {   //multipliers from first area
                if(j!=p) {  //without singularity
                    delta=gsl_complex_sqrt_real(-pow(k,2)+pow(MMN[j]/a,2));
                    if(GSL_REAL(delta)==0) delta=gsl_complex_mul_real(delta,-1);
                    z=gsl_complex_add(beta1,beta);
                    z=gsl_complex_exp(gsl_complex_mul_real(z,-a/M_PI/(j+1)));
                    stream << str1.setNum(j+1) << "\t" << str.setNum(gsl_complex_abs(z)) << "\t";
                    z=gsl_complex_mul(z,gsl_complex_div(gsl_complex_add(beta1,delta),gsl_complex_sub(beta,delta)));
                    Aq=gsl_complex_mul(Aq,z);
                    stream << str.setNum(gsl_complex_abs(Aq)) << "\n";
                }
                //! superfluity
                delta=gsl_complex_sqrt_real(-pow(k,2)+pow(MMN[itM]/Shield,2));
                if(GSL_REAL(delta)==0) sign=-1; else sign=1;
                //delta=gsl_complex_sqrt(gsl_complex_add_real(eps,pow(MMN[itM]/Shield,2)));
                delta=gsl_complex_mul_real(delta,sign);/*if(GSL_REAL(delta)==0)*/
                z=gsl_complex_add(beta1,beta);
                z=gsl_complex_exp(gsl_complex_mul_real(z,Shield/M_PI/(itM+1)));
                stream << str1.setNum(itM+1) << "\t" << str.setNum(gsl_complex_abs(z)) << "\t";
                z=gsl_complex_mul(z,gsl_complex_div(gsl_complex_sub(beta,delta),gsl_complex_add(beta1,delta)));
                Aq=gsl_complex_mul(Aq,z);
                stream << str.setNum(gsl_complex_abs(Aq)) << "\n";
                itM++;
                for(int m=0; m<difQ[3]-1; m++)
                {
                    //!  coax
                    delta=gsl_complex_sqrt_real(-pow(k,2)+pow(HMN[itH]/a,2));
                    if(GSL_REAL(delta)==0) sign=-1; else sign=1;
                    //delta=gsl_complex_sqrt(gsl_complex_add_real(eps,pow(HMN[itH]/a,2)));
                    delta=gsl_complex_mul_real(delta,sign); /*if(GSL_REAL(delta)==0)*/
                    //if (p==0) stack->append("Ro:\t"+str.setNum(GSL_REAL(delta))+"\t"+str1.setNum(GSL_IMAG(delta)));
                    z=gsl_complex_add(beta1,beta);
                    z=gsl_complex_exp(gsl_complex_mul_real(z,-(Shield-a)/M_PI/(itH+1)));
                    stream << str1.setNum(itH+1) << "\t" << str.setNum(gsl_complex_abs(z)) << "\t";
                    z=gsl_complex_mul(z,gsl_complex_div(gsl_complex_add(beta1,delta),gsl_complex_sub(beta,delta)));
                    Aq=gsl_complex_mul(Aq,z);
                    stream << str.setNum(gsl_complex_abs(Aq)) << "\n";
                    itH++;
                    //! large circular WG
                    delta=gsl_complex_sqrt_real(-pow(k,2)+pow(MMN[itM]/Shield,2));
                    if(GSL_REAL(delta)==0) sign=-1; else sign=1;
                    //delta=gsl_complex_sqrt(gsl_complex_add_real(eps,pow(MMN[itM]/Shield,2)));
                    delta=gsl_complex_mul_real(delta,sign);/*if(GSL_REAL(delta)==0)*/
                    z=gsl_complex_add(beta1,beta);
                    z=gsl_complex_exp(gsl_complex_mul_real(z,Shield/M_PI/(itM+1)));
                    stream << str1.setNum(itM+1) << "\t" << str.setNum(gsl_complex_abs(z)) << "\t";
                    z=gsl_complex_mul(z,gsl_complex_div(gsl_complex_sub(beta,delta),gsl_complex_add(beta1,delta)));
                    Aq=gsl_complex_mul(Aq,z);
                    stream << str.setNum(gsl_complex_abs(Aq)) << "\n";
                    itM++;
                }
            }
            stack->append("first\t"+str.setNum(gsl_complex_abs(Aq))+"\t"+str1.setNum(gsl_complex_arg(Aq)));
            stream << "A[" << str.setNum(p+1) << "]:\n";
            stream << str.setNum(GSL_REAL(Aq)) << "\t" << str1.setNum(GSL_IMAG(Aq)) << "\n\n";
            gsl_vector_complex_set(Xnew,p,Aq);
        }
        Refract=gsl_vector_complex_get(Xnew,0);
        writefile.close();
    }
}

void Vainshtein::analiticE(QTextBrowser *stack)
{
    int N,sign=1;
    QString str,str1;
    QFile writefile("ANALITIC_E.CSV");
    QTextStream stream(&writefile);
    //initialisation
    difQ[2]=difQ[1]*(difQ[3]-1);//! convergence
    N=difQ[1]+difQ[2];    //without TEM
    for (int j=80; j<=difQ[2]; j++) EMN[j-1]=zerosCoaxE(mI,j);
    for (int j=80; j<=N; j++) VMN[j-1]=gsl_sf_bessel_zero_Jnu(mI,j);
    gsl_complex z,Aq,Aq2,beta1,beta,delta;
    gsl_vector_complex * Xnew = gsl_vector_complex_alloc (N);
    //eps=gsl_complex_pow_real(gsl_complex_mul_real(eps,-k),2);
    eps=gsl_complex_mul_real(eps,pow(k,2));
    beta1=gsl_complex_sqrt_real(-pow(k,2)+pow(VMN[0]/a,2));
    beta1=gsl_complex_mul_real(beta1,-1);
    if(writefile.open(QIODevice::WriteOnly))
    {
        for(int p=0; p<difQ[1]; p++)
        {
            //GSL_SET_COMPLEX(&beta,1,p+1);
            beta=gsl_complex_sqrt_real(-pow(k,2)+pow(VMN[p]/a,2));
            //GSL_SET_COMPLEX(&Aq,gsl_sf_bessel_Jn(2,VMN[0])/gsl_sf_bessel_Jn(2,VMN[p]),0);
            if(GSL_REAL(beta)==0) beta=gsl_complex_mul_real(beta,-1);
            GSL_SET_COMPLEX(&z,((a*log(a/(Shield-a))-Shield*log(Shield/(Shield-a)))-a/(p+1))/M_PI,0);
            //GSL_SET_COMPLEX(&z,((a*log((double)difQ[1]/difQ[2])+Shield*log((double)difQ[2]/(difQ[2]+difQ[1])))-a/(p+1))/M_PI,0);
            z=gsl_complex_mul(z,gsl_complex_add(beta1,beta));
            Aq=gsl_complex_mul_real(gsl_complex_exp(z),gsl_sf_bessel_Jn(2,VMN[0])/gsl_sf_bessel_Jn(2,VMN[p]));
            GSL_SET_COMPLEX(&Aq2,gsl_sf_bessel_Jn(2,VMN[0])/gsl_sf_bessel_Jn(2,VMN[p]),0);
            stream << "p \t exp \t div \t Aq \t Aq2 \n";
            for(int j=0; j<N; j++)
            {
                delta=gsl_complex_sqrt_real(-pow(k,2)+pow(VMN[j]/Shield,2));
                if(GSL_REAL(delta)==0) sign=-1; else sign=1;
                //delta=gsl_complex_sqrt(gsl_complex_sub_real(eps,pow(VMN[j]/Shield,2)));
                /*if(GSL_REAL(delta)==0)*/ delta=gsl_complex_mul_real(delta,sign);
                z=gsl_complex_add(beta1,beta);
                z=gsl_complex_exp(gsl_complex_mul_real(z,Shield/M_PI/(j+1)));
                stream << str1.setNum(j+1) << "\t" << str.setNum(gsl_complex_abs(z)) << "\t";
                z=gsl_complex_mul(z,gsl_complex_div(gsl_complex_sub(beta,delta),gsl_complex_add(beta1,delta)));
                Aq=gsl_complex_mul(Aq,z);
                z=gsl_complex_div(gsl_complex_sub(beta,delta),gsl_complex_add(beta1,delta));
                stream << str.setNum(gsl_complex_abs(z)) << "\t";
                Aq2=gsl_complex_mul(Aq2,z);
                stream << str.setNum(gsl_complex_abs(Aq)) << "\t";
                stream << str.setNum(gsl_complex_abs(Aq2)) << "\n";
            }
            stream << "n \t exp \t div \t Aq \t Aq2 \n";
            for(int j=0; j<difQ[1]; j++)
            {
                if(j!=p) {
                    delta=gsl_complex_sqrt_real(-pow(k,2)+pow(VMN[j]/a,2));
                    if(GSL_REAL(delta)==0) delta=gsl_complex_mul_real(delta,-1);
                    z=gsl_complex_add(beta1,beta);
                    z=gsl_complex_exp(gsl_complex_mul_real(z,-a/M_PI/(j+1)));
                    stream << str1.setNum(j+1) << "\t" << str.setNum(gsl_complex_abs(z)) << "\t";
                    z=gsl_complex_mul(z,gsl_complex_div(gsl_complex_add(beta1,delta),gsl_complex_sub(beta,delta)));
                    Aq=gsl_complex_mul(Aq,z);
                    z=gsl_complex_div(gsl_complex_add(beta1,delta),gsl_complex_sub(beta,delta));
                    stream << str.setNum(gsl_complex_abs(z)) << "\t";
                    Aq2=gsl_complex_mul(Aq2,z);
                    stream << str.setNum(gsl_complex_abs(Aq)) << "\t";
                    stream << str.setNum(gsl_complex_abs(Aq2)) << "\n";
                } }
            stream << "m \t exp \t div \t Aq \t Aq2 \n";
            for(int j=0; j<difQ[2]; j++)
            {
                delta=gsl_complex_sqrt_real(-pow(k,2)+pow(EMN[j]/a,2));
                if(GSL_REAL(delta)==0) sign=-1; else sign=1;
                //delta=gsl_complex_sqrt(gsl_complex_sub_real(eps,pow(EMN[j]/a,2)));
                /*if(GSL_REAL(delta)==0)*/ delta=gsl_complex_mul_real(delta,sign);
                //if (p==0) stack->append("Ro:\t"+str.setNum(GSL_REAL(delta))+"\t"+str1.setNum(GSL_IMAG(delta)));
                z=gsl_complex_add(beta1,beta);
                z=gsl_complex_exp(gsl_complex_mul_real(z,-(Shield-a)/M_PI/(j+1)));
                stream << str1.setNum(j+1) << "\t" << str.setNum(gsl_complex_abs(z)) << "\t";
                z=gsl_complex_mul(z,gsl_complex_div(gsl_complex_add(beta1,delta),gsl_complex_sub(beta,delta)));
                Aq=gsl_complex_mul(Aq,z);
                z=gsl_complex_div(gsl_complex_add(beta1,delta),gsl_complex_sub(beta,delta));
                stream << str.setNum(gsl_complex_abs(z)) << "\t";
                Aq2=gsl_complex_mul(Aq2,z);
                stream << str.setNum(gsl_complex_abs(Aq)) << "\t";
                stream << str.setNum(gsl_complex_abs(Aq2)) << "\n";
            }
            //stack->append(str.setNum(GSL_REAL(Aq))+"+i"+str1.setNum(GSL_IMAG(Aq)));
            stack->append("first\t"+str.setNum(gsl_complex_abs(Aq))+"\t"+str1.setNum(gsl_complex_arg(Aq)));
            stack->append("second\t"+str.setNum(gsl_complex_abs(Aq2))+"\t"+str1.setNum(gsl_complex_arg(Aq2)));
            stream << "A[" << str.setNum(p+1) << "]:\n";
            stream << str.setNum(GSL_REAL(Aq)) << "\t" << str1.setNum(GSL_IMAG(Aq)) << "\t";
            stream << str.setNum(GSL_REAL(Aq2)) << "\t" << str1.setNum(GSL_IMAG(Aq2)) << "\n\n";
            gsl_vector_complex_set(Xnew,p,Aq);
        }
        Refract=gsl_vector_complex_get(Xnew,0);
        writefile.close();
    }
}


void Vainshtein::withoutInt2(QTextBrowser *stack)
{
    int N;
    QString str,str1;
    QFile writefile("SOLVE.CSV");
    QTextStream stream(&writefile);
    //initialisation
    N=difQ[1]+difQ[2];    //without TEM
    gsl_complex z;
    gsl_matrix_complex * A = gsl_matrix_complex_alloc (N, N);
    gsl_vector_complex * B = gsl_vector_complex_alloc (N);
    gsl_vector_complex * X = gsl_vector_complex_alloc (N);
    int s;
    gsl_permutation * P = gsl_permutation_alloc (N);
    if(writefile.open(QIODevice::WriteOnly))
    {
        stack->append("k="+str.setNum(k));
        stack->append("N="+str.setNum(N));
        //stack->append("J2="+str.setNum(gsl_sf_bessel_Jn(2,MMN[0])));
        //end of initialisation
        gsl_matrix_complex_set_zero (A);
        GSL_SET_COMPLEX(&z,1,0);
        gsl_vector_complex_set_all (B,z);
        //! \Fn coefficients
        gsl_complex beta,beta1,delta;
        eps=gsl_complex_pow_real(gsl_complex_mul_real(eps,k),2);
        beta1=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[0]/a,2));
        double Ji;
        for (int j=0; j<N; j++)	//! first area with R1=1
        {
            delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/Shield,2));
            //delta=gsl_complex_sqrt(gsl_complex_sub_real(eps,pow(MMN[j]/Shield,2)));
            for (int q=0; q<difQ[1]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/a,2));
                Ji=gsl_sf_bessel_Jn(2,MMN[q])/gsl_sf_bessel_Jn(2,MMN[0]);
                z=gsl_complex_add(beta1,delta);
                z=gsl_complex_div(z,gsl_complex_sub(beta,delta));
                gsl_matrix_complex_set (A,j,q,gsl_complex_mul_real(z,Ji));
            }
            for (int q=0; q<difQ[2]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(HMN[q]/a,2));
                //beta=gsl_complex_sqrt(gsl_complex_sub_real(eps,pow(HMN[q]/a,2)));
                Ji=gsl_sf_bessel_Jn(1,HMN[q])*gsl_sf_bessel_Yn(2,HMN[q])-gsl_sf_bessel_Yn(1,HMN[q])*gsl_sf_bessel_Jn(2,HMN[q]);
                Ji/=gsl_sf_bessel_Jn(2,MMN[0]);
                z=gsl_complex_add(beta1,delta);
                z=gsl_complex_div(z,gsl_complex_sub(beta,delta));
                gsl_matrix_complex_set (A,j,difQ[1]+q,gsl_complex_mul_real(z,Ji));
            }
            /*! Without norming to "1" in left part  */ /*
//            gsl_vector_complex_set (B,j,gsl_complex_mul_real(z,Ji*MMN[0]/a));
//            delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/Shield,2));
//            alph=pow(MMN[0]/a,2)-pow(MMN[j]/Shield,2);
//            for (int q=0; q<difQ[1]; q++)
//            {
//                beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/a,2));
//                Ji=gsl_sf_bessel_Jn(2,MMN[q])/(pow(MMN[q]/a,2)-pow(MMN[j]/Shield,2))*alph;
//                z=gsl_complex_add(beta,delta);
//                gsl_matrix_complex_set (A,j,q,gsl_complex_mul_real(z,Ji*MMN[q]/a));
//            }
//            for (int q=0; q<difQ[2]; q++)
//            {
//                beta=gsl_complex_sqrt_real(pow(k,2)-pow(HMN[q]/a,2));
//                Ji=gsl_sf_bessel_Jn(1,HMN[q])*gsl_sf_bessel_Yn(2,HMN[q])-gsl_sf_bessel_Yn(1,HMN[q])*gsl_sf_bessel_Jn(2,HMN[q]);
//                Ji*=alph/(pow(HMN[q]/a,2)-pow(MMN[j]/Shield,2));
//                z=gsl_complex_add(beta,delta);
//                gsl_matrix_complex_set (A,j,difQ[1]+q,gsl_complex_mul_real(z,Ji*HMN[q]/a));
//            }
//            beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[0]/a,2));
//            Ji=gsl_sf_bessel_Jn(2,MMN[0]);
//            z=gsl_complex_sub(beta,delta);
//            gsl_vector_complex_set (B,j,gsl_complex_mul_real(z,Ji*MMN[0]/a));
              */
        }
        //! \Fn solving
        gsl_matrix_complex * Acopy = gsl_matrix_complex_alloc (N, N);
        gsl_matrix_complex_memcpy(Acopy,A);
        gsl_linalg_complex_LU_decomp (A, P, &s);
        gsl_linalg_complex_LU_solve (A, P, B, X);
        viewinfo (&stream, N, Acopy, B, X, P);
        //! checking
        stream << "Zeros MMN: \n";
        gsl_vector_complex * Check = gsl_vector_complex_alloc(N);
        for (int j=0; j<N; j++)
        {
            stream << str.setNum(MMN[j]) << endl;
        }
        stream << "Checking: \n";
        gsl_vector_complex_set_zero(Check);
        for (int j=0; j<N; j++)
        {
            GSL_SET_COMPLEX(&z,0,0);
            for (int q=0; q<N; q++)
            {
                z=gsl_complex_add(z,gsl_complex_mul(gsl_matrix_complex_get(Acopy,j,q),gsl_vector_complex_get(X,q)));
            }
            gsl_vector_complex_set(Check,j,z);
            stream << str.setNum(GSL_REAL(gsl_vector_complex_get(Check,j))) << "\t" << str1.setNum(GSL_IMAG(gsl_vector_complex_get(Check,j))) << endl;
        }
    }
    writefile.close();
    for(int j=0;j<difQ[1];j++)
    {
        z=gsl_vector_complex_get(X,j);
        str=str.setNum(GSL_REAL(z));
        str1=str1.setNum(GSL_IMAG(z));
        //stack->append(str+" + i "+str1);
        stack->append("Coef: "+str.setNum(gsl_complex_abs(z))+"\t"+str1.setNum(gsl_complex_arg(z)));
    }
    Refract=gsl_vector_complex_get(X,0);
    gsl_permutation_free (P);
}

void Vainshtein::withoutInt(QTextBrowser *stack)
{
    int q,N;
    QString str,str1;
    QFile writefile("SOLVE.CSV");
    QTextStream stream(&writefile);
    //initialisation
    //N=difQ[1]+difQ[2];    //without TEM
    N=difQ[1];    //without TEM
    gsl_complex z;
    gsl_matrix_complex * A = gsl_matrix_complex_alloc (N, N);
    gsl_vector_complex * B = gsl_vector_complex_alloc (N);
    gsl_vector_complex * X = gsl_vector_complex_alloc (N);
    int s;
    gsl_permutation * P = gsl_permutation_alloc (N);
    if(writefile.open(QIODevice::WriteOnly))
    {
        stack->append("k="+str.setNum(k));
        stack->append("N="+str.setNum(N));
        //end of initialisation
        gsl_matrix_complex_set_zero (A);
        gsl_vector_complex_set_zero (B);
        GSL_SET_COMPLEX(&z,1,1);
        //! \Fn coefficients
        gsl_complex beta,delta;
        double Ji;
        for (int j=0; j<N; j++)	//! first area with R1=1
        {
            delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/Shield,2));
            Ji=gsl_sf_bessel_Jn(1,MMN[j]/Shield*a);
            for (q=0; q<difQ[1]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/a,2));
                z=gsl_complex_add(beta,delta);
                gsl_matrix_complex_set (A,q,j,gsl_complex_mul_real(z,Ji/(pow(MMN[q]/a,2)-pow(MMN[j]/Shield,2))));
            }
        }
        beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[0]/a,2));
        Ji=pow(besselS(1,MMN[0])*a,2)/gsl_sf_bessel_Jn(2,MMN[0]);
        gsl_vector_complex_set (B,0,gsl_complex_mul_real(beta,Ji/MMN[0]));
        //! \Fn solving
        gsl_matrix_complex * Acopy = gsl_matrix_complex_alloc (N, N);
        gsl_matrix_complex_memcpy(Acopy,A);
        gsl_linalg_complex_LU_decomp (A, P, &s);
        gsl_linalg_complex_LU_solve (A, P, B, X);
        viewinfo (&stream, N, Acopy, B, X, P);
        gsl_complex coefA;  //itog
        stream << "vector AA:\n";
        for (int q=0; q<difQ[1]; q++)
        {
            beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/a,2));
            for (int j=0; j<N; j++)
            {
                delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/Shield,2));
                z=gsl_complex_div(gsl_vector_complex_get (X,j),gsl_complex_sub(beta,delta));
                z=gsl_complex_mul_real(coefA,gsl_sf_bessel_Jn(1,MMN[j]/Shield*a));
                coefA=gsl_complex_add(coefA,z);
            }
            coefA=gsl_complex_div(coefA,beta);
            coefA=gsl_complex_div_real(coefA,pow(besselS(1,MMN[q])*a,2)/MMN[q]/gsl_sf_bessel_Jn(2,MMN[q]));
            stream << str.setNum(GSL_REAL(coefA)) << "\t" << str1.setNum(GSL_IMAG(coefA)) << endl;
            stack->append("Coef "+str1.setNum(q)+": "+str.setNum(gsl_complex_abs(coefA)));
        }
        //! checking
        stream << "Zeros MMN: \n";
        gsl_vector_complex * Check = gsl_vector_complex_alloc(N);
        for (int j=0; j<N; j++)
        {
            stream << str.setNum(MMN[j]) << endl;
        }
        stream << "Checking: \n";
        gsl_vector_complex_set_zero(Check);
        for (int j=0; j<N; j++)
        {
            GSL_SET_COMPLEX(&z,0,0);
            for (int q=0; q<N; q++)
            {
                z=gsl_complex_add(z,gsl_complex_mul(gsl_matrix_complex_get(Acopy,j,q),gsl_vector_complex_get(X,q)));
            }
            gsl_vector_complex_set(Check,j,z);
            stream << str.setNum(GSL_REAL(gsl_vector_complex_get(Check,j))) << "\t" << str1.setNum(GSL_IMAG(gsl_vector_complex_get(Check,j))) << endl;
        }
    }
    writefile.close();
    gsl_permutation_free (P);
}

//! \fn Diffraction issue
void Vainshtein::diffraction(QTextBrowser *stack)
{
    int N;
    QString str,str1;
    QFile writefile("SOLVE.CSV");
    QTextStream stream(&writefile);
    //initialisation
    N=difQ[1]+difQ[2];    //without TEM
    gsl_complex z;
    gsl_matrix_complex * A = gsl_matrix_complex_alloc (N, N);
    gsl_vector_complex * B = gsl_vector_complex_alloc (N);
    gsl_vector_complex * X = gsl_vector_complex_alloc (N);
    int s;
    gsl_permutation * P = gsl_permutation_alloc (N);
    if(writefile.open(QIODevice::WriteOnly))
    {
        stack->append("k="+str.setNum(k));
        stack->append("N="+str.setNum(N));
        //end of initialisation
        gsl_matrix_complex_set_zero (A);
        GSL_SET_COMPLEX(&z,1,1);
        //! \Fn coefficients
        gsl_complex beta,delta;
        double Ji;
        stream << "first area" << endl;
        for (int j=0; j<N; j++)	//! first area with R1=1
        {
            delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/Shield,2));
            stream << "delta=\t" << str.setNum(GSL_REAL(delta)) << " + i*" << str1.setNum(GSL_IMAG(delta)) << endl;
            stream << "alpha[" << str1.setNum(j) << "] = \t" << str.setNum(MMN[j]/a) << endl;
            for (int q=0; q<difQ[1]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/a,2));
                Ji=amq(MMN[q],MMN[j]);
                stream << "Ji[" << str1.setNum(q) << "] = \t" << str.setNum(Ji) << endl;
                z=gsl_complex_add(beta,delta);
                stream << "Sbeta=\t" << str.setNum(GSL_REAL(z)) << " + i*" << str1.setNum(GSL_IMAG(z)) << endl;
                gsl_matrix_complex_set (A,j,q,gsl_complex_mul_real(z,Ji));
            }
            for (int q=0; q<difQ[2]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(HMN[q],2));
                Ji=bnq(HMN[q],MMN[j]);
                stream << "Ji[" << str1.setNum(q) << "] = \t" << str.setNum(Ji) << endl;
                z=gsl_complex_add(beta,delta);
                stream << "Sbeta=\t" << str.setNum(GSL_REAL(z)) << " + i*" << str1.setNum(GSL_IMAG(z)) << endl;
                gsl_matrix_complex_set (A,j,difQ[1]+q,gsl_complex_mul_real(z,Ji));
            }
            beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[nI-1]/a,2));
            Ji=amq(MMN[nI-1],MMN[j]);
            stream << "Ji[" << str1.setNum(j) << "] = \t" << str.setNum(Ji) << endl;
            z=gsl_complex_sub(beta,delta);
            stream << "Rbeta=\t" << str.setNum(GSL_REAL(z)) << " + i*" << str1.setNum(GSL_IMAG(z)) << endl;
            gsl_vector_complex_set (B,j,gsl_complex_mul_real(z,Ji));
        }
        //! \Fn solving
        gsl_matrix_complex * Acopy = gsl_matrix_complex_alloc (N, N);
        gsl_matrix_complex_memcpy(Acopy,A);
        gsl_linalg_complex_LU_decomp (A, P, &s);
        gsl_linalg_complex_LU_solve (A, P, B, X);
        viewinfo (&stream, N, Acopy, B, X, P);
        //! checking
        stream << "Zeros MMN: \n";
        gsl_vector_complex * Check = gsl_vector_complex_alloc(N);
        for (int j=0; j<N; j++)
        {
            stream << str.setNum(MMN[j]) << endl;
        }
        stream << "Checking: \n";
        gsl_vector_complex_set_zero(Check);
        for (int j=0; j<N; j++)
        {
            GSL_SET_COMPLEX(&z,0,0);
            for (int q=0; q<N; q++)
            {
                z=gsl_complex_add(z,gsl_complex_mul(gsl_matrix_complex_get(Acopy,j,q),gsl_vector_complex_get(X,q)));
            }
            gsl_vector_complex_set(Check,j,z);
            stream << str.setNum(GSL_REAL(gsl_vector_complex_get(Check,j))) << "\t" << str1.setNum(GSL_IMAG(gsl_vector_complex_get(Check,j))) << endl;
        }
    }
    writefile.close();
    for(int j=0;j<N;j++)
    {
        z=gsl_vector_complex_get(X,j);
        str=str.setNum(GSL_REAL(z));
        str1=str1.setNum(GSL_IMAG(z));
        stack->append(str+" + i "+str1);
        stack->append("Coef "+str1.setNum(j)+": "+str.setNum(gsl_complex_abs(z)));
    }
    gsl_permutation_free (P);
}

void Vainshtein::diff_new(QTextBrowser *stack)
{//! new method
    int q,N;
    QString str,str1;
    QFile writefile("SOLVE.CSV");
    QTextStream stream(&writefile);
    //initialisation
    N=difQ[1]+difQ[2];    //without TEM
    //N=difQ[1];    //without TEM
    gsl_complex z;
    gsl_matrix_complex * A = gsl_matrix_complex_alloc (N, N);
    gsl_vector_complex * B = gsl_vector_complex_alloc (N);
    gsl_vector_complex * X = gsl_vector_complex_alloc (N);
    int s;
    gsl_permutation * P = gsl_permutation_alloc (N);
    if(writefile.open(QIODevice::WriteOnly))
    {
        stack->append("k="+str.setNum(k));
        stack->append("N="+str.setNum(N));
        //end of initialisation
        gsl_matrix_complex_set_zero (A);
        gsl_vector_complex_set_zero (B);
        GSL_SET_COMPLEX(&z,1,1);
        //! \Fn coefficients
        gsl_complex beta,delta;
        double Ji;
        stream << "first area" << endl;
        for (int j=0; j<N; j++)	//! first area with R1=1
        {
            delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/Shield,2));
            stream << "delta=\t" << str.setNum(GSL_REAL(delta)) << " + i*" << str1.setNum(GSL_IMAG(delta)) << endl;
            stream << "alpha[" << str1.setNum(j) << "] = \t" << str.setNum(MMN[j]/a) << endl;
            for (q=0; q<difQ[1]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/a,2));
                Ji=amq(MMN[q],MMN[j]);
                stream << "Ji[" << str1.setNum(q) << "] = \t" << str.setNum(Ji) << endl;
                z=gsl_complex_add(beta,delta);
                stream << "Sbeta=\t" << str.setNum(GSL_REAL(z)) << " + i*" << str1.setNum(GSL_IMAG(z)) << endl;
                gsl_matrix_complex_set (A,q,j,gsl_complex_mul_real(z,Ji));
            }
            for (q=0; q<difQ[2]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(HMN[q],2));
                Ji=bnq(HMN[q],MMN[j]);
                stream << "Ji[" << str1.setNum(q) << "] = \t" << str.setNum(Ji) << endl;
                z=gsl_complex_add(beta,delta);
                stream << "Sbeta=\t" << str.setNum(GSL_REAL(z)) << " + i*" << str1.setNum(GSL_IMAG(z)) << endl;
                gsl_matrix_complex_set (A,difQ[1]+q,j,gsl_complex_mul_real(z,Ji));
            }
        }
        beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[0]/a,2));
        Ji=pow(besselS(1,MMN[0])*a,2);
        stream << "Ji = \t" << str.setNum(Ji) << endl;
        stream << "Rbeta=\t" << str.setNum(GSL_REAL(beta)) << " + i*" << str1.setNum(GSL_IMAG(beta)) << endl;
        gsl_vector_complex_set (B,0,gsl_complex_mul_real(beta,Ji));
        //! \Fn solving
        gsl_matrix_complex * Acopy = gsl_matrix_complex_alloc (N, N);
        gsl_matrix_complex_memcpy(Acopy,A);
        gsl_linalg_complex_LU_decomp (A, P, &s);
        gsl_linalg_complex_LU_solve (A, P, B, X);
        viewinfo (&stream, N, Acopy, B, X, P);
        gsl_complex coefA;  //itog
        stream << "vector AA:\n";
        for (int q=0; q<difQ[1]; q++)
        {
            beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/a,2));
            for (int j=0; j<N; j++)
            {
                delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/Shield,2));
                Ji=amq(MMN[q],MMN[j]);
                z=gsl_complex_sub(beta,delta);
                z=gsl_complex_mul_real(z,Ji);
                z=gsl_complex_mul(z,gsl_vector_complex_get (X,j));
                coefA=gsl_complex_add(coefA,z);
            }
            coefA=gsl_complex_div(coefA,beta);
            coefA=gsl_complex_div_real(coefA,pow(besselS(1,MMN[q])*a,2));
            stream << str.setNum(GSL_REAL(coefA)) << "\t" << str1.setNum(GSL_IMAG(coefA)) << endl;
            stack->append("Coef "+str1.setNum(q)+": "+str.setNum(gsl_complex_abs(coefA)));
        }
        //! checking
        stream << "Zeros MMN: \n";
        gsl_vector_complex * Check = gsl_vector_complex_alloc(N);
        for (int j=0; j<N; j++)
        {
            stream << str.setNum(MMN[j]) << endl;
        }
        stream << "Checking: \n";
        gsl_vector_complex_set_zero(Check);
        for (int j=0; j<N; j++)
        {
            GSL_SET_COMPLEX(&z,0,0);
            for (int q=0; q<N; q++)
            {
                z=gsl_complex_add(z,gsl_complex_mul(gsl_matrix_complex_get(Acopy,j,q),gsl_vector_complex_get(X,q)));
            }
            gsl_vector_complex_set(Check,j,z);
            stream << str.setNum(GSL_REAL(gsl_vector_complex_get(Check,j))) << "\t" << str1.setNum(GSL_IMAG(gsl_vector_complex_get(Check,j))) << endl;
        }
    }
    writefile.close();
    gsl_permutation_free (P);
}

void Vainshtein::wenger(QTextBrowser *stack)
{
    int N;
    QString str,str1;
    //QFile writefile("SOLVE.CSV");
    //QTextStream stream(&writefile);
    N=difQ[1]+difQ[2];    //without TEM
    Shield=a;
    a=Shield/1.5;
    for (int j=1; j<=40; j++) HMN[j-1]=zerosCoaxH(mI,j);
    gsl_complex z,Res,delta;
    double beta;
    stack->append("k="+str.setNum(k));
    stack->append("N="+str.setNum(N));
    stack->append("II:\t"+str.setNum(pow(k,2)-pow(MMN[0]/a,2)));
    stack->append("III:\t"+str.setNum(pow(k,2)-pow(HMN[0]/a,2)));
    beta=sqrt(pow(k,2)-pow(MMN[0]/Shield,2));
    GSL_SET_COMPLEX(&z,0,-beta/M_PI*2*(a*log(a/(Shield-a))-Shield*log(Shield/(Shield-a))));
    Res=gsl_complex_exp(z);
    Res=gsl_complex_mul_real(Res,-1);
    for (int j=1; j<difQ[1]; j++)
    {
        GSL_SET_COMPLEX(&delta,sqrt(pow(MMN[j]/Shield,2)-pow(k,2)),beta);
        Res=gsl_complex_mul(Res,gsl_complex_div_real(gsl_complex_pow_real(delta,2),pow(gsl_complex_abs(delta),2)));
    }
    for (int j=0; j<difQ[1]; j++)
    {
        GSL_SET_COMPLEX(&delta,sqrt(pow(MMN[j]/a,2)-pow(k,2)),-beta);
        Res=gsl_complex_mul(Res,gsl_complex_div_real(gsl_complex_pow_real(delta,2),pow(gsl_complex_abs(delta),2)));
    }
    for (int j=0; j<difQ[2]; j++)
    {
        GSL_SET_COMPLEX(&delta,sqrt(pow(HMN[j],2)-pow(k,2)),-beta);
        Res=gsl_complex_mul(Res,gsl_complex_div_real(gsl_complex_pow_real(delta,2),pow(gsl_complex_abs(delta),2)));
    }
    stack->append("Coef=\t"+str1.setNum(GSL_REAL(Res))+"+i"+str.setNum(GSL_IMAG(Res)));
    stack->append("|Coef|:\t"+str.setNum(gsl_complex_abs(Res)));
    stack->append("arg(Coef):\t"+str.setNum(gsl_complex_arg(Res)));
}

void Vainshtein::diff_check(QTextBrowser *stack)
{//! new method
    int q,N;
    QString str,str1;
    QFile writefile("SOLVE.CSV");
    QTextStream stream(&writefile);
    //initialisation
    N=difQ[1]+difQ[2];    //without TEM
    Shield=a;
    a=Shield/1.75;
    for (int j=1; j<=10; j++) HMN[j-1]=zerosCoaxH(mI,j);
    //N=difQ[1];    //without TEM
    gsl_complex z;
    gsl_matrix_complex * A = gsl_matrix_complex_alloc (N, N);
    gsl_vector_complex * B = gsl_vector_complex_alloc (N);
    gsl_vector_complex * X = gsl_vector_complex_alloc (N);
    int s;
    gsl_permutation * P = gsl_permutation_alloc (N);
    if(writefile.open(QIODevice::WriteOnly))
    {
        stack->append("k="+str.setNum(k));
        stack->append("N="+str.setNum(N));
        stack->append("II:\t"+str.setNum(pow(k,2)-pow(MMN[0]/a,2)));
        stack->append("III:\t"+str.setNum(pow(k,2)-pow(HMN[0]/a,2)));
        //end of initialisation
        gsl_matrix_complex_set_zero (A);
        gsl_vector_complex_set_zero (B);
        GSL_SET_COMPLEX(&z,1,1);
        //! \Fn coefficients
        gsl_complex beta,delta;
        double Ji;
        stream << "left part" << endl;
        for (int j=0; j<N; j++)	//! first area with R1=1
        {
            delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/Shield,2));
            stream << "delta=\t" << str.setNum(GSL_REAL(delta)) << " + i*" << str1.setNum(GSL_IMAG(delta)) << endl;
            stream << "alpha[" << str1.setNum(j) << "] = \t" << str.setNum(MMN[j]/a) << endl;
            for (q=0; q<difQ[1]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/a,2));
                Ji=amq(MMN[q],MMN[j]);
                stream << "Ji[" << str1.setNum(q) << "] = \t" << str.setNum(Ji) << endl;
                z=gsl_complex_add(beta,delta);
                stream << "Sbeta=\t" << str.setNum(GSL_REAL(z)) << " + i*" << str1.setNum(GSL_IMAG(z)) << endl;
                gsl_matrix_complex_set (A,j,q,gsl_complex_mul_real(z,Ji));
            }
            for (q=0; q<difQ[2]; q++)
            {
                beta=gsl_complex_sqrt_real(pow(k,2)-pow(HMN[q],2));
                Ji=bnq(HMN[q],MMN[j]);
                stream << "Ji[" << str1.setNum(q) << "] = \t" << str.setNum(Ji) << endl;
                z=gsl_complex_add(beta,delta);
                stream << "Sbeta=\t" << str.setNum(GSL_REAL(z)) << " + i*" << str1.setNum(GSL_IMAG(z)) << endl;
                gsl_matrix_complex_set (A,j,difQ[1]+q,gsl_complex_mul_real(z,Ji));
            }
        }
        beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[0]/Shield,2));
        Ji=pow(besselS(1,MMN[0])*Shield,2);
        stream << "Ji = \t" << str.setNum(Ji) << endl;
        stream << "Rbeta=\t" << str.setNum(GSL_REAL(beta)) << " + i*" << str1.setNum(GSL_IMAG(beta)) << endl;
        gsl_vector_complex_set (B,0,gsl_complex_mul_real(beta,Ji));
        //! \Fn solving
        gsl_matrix_complex * Acopy = gsl_matrix_complex_alloc (N, N);
        gsl_matrix_complex_memcpy(Acopy,A);
        gsl_linalg_complex_LU_decomp (A, P, &s);
        gsl_linalg_complex_LU_solve (A, P, B, X);
        viewinfo (&stream, N, Acopy, B, X, P);
        gsl_complex coefA;  //itog
        stream << "vector AA:\n";
        for (int q=0; q<N; q++)
        {
            beta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[q]/Shield,2));
            for (int j=0; j<difQ[1]; j++)
            {
                delta=gsl_complex_sqrt_real(pow(k,2)-pow(MMN[j]/a,2));
                Ji=amq(MMN[j],MMN[q]);
                z=gsl_complex_sub(beta,delta);
                z=gsl_complex_mul_real(z,Ji);
                z=gsl_complex_mul(z,gsl_vector_complex_get (X,j));
                coefA=gsl_complex_add(coefA,z);
            }
            for (int j=0; j<difQ[2]; j++)
            {
                delta=gsl_complex_sqrt_real(pow(k,2)-pow(HMN[j],2));
                Ji=bnq(HMN[j],MMN[q]);
                z=gsl_complex_sub(beta,delta);
                z=gsl_complex_mul_real(z,Ji);
                z=gsl_complex_mul(z,gsl_vector_complex_get (X,j));
                coefA=gsl_complex_add(coefA,z);
            }
            //coefA=gsl_complex_div(coefA,beta);
            //coefA=gsl_complex_div_real(coefA,pow(besselS(1,MMN[q])*Shield,2));
            stream << str.setNum(GSL_REAL(coefA)) << "\t" << str1.setNum(GSL_IMAG(coefA)) << endl;
            stack->append("Coef "+str1.setNum(q)+": "+str.setNum(gsl_complex_abs(coefA)));
        }
        //! checking
        stream << "Zeros HMN: \n";
        gsl_vector_complex * Check = gsl_vector_complex_alloc(N);
        for (int j=0; j<difQ[2]; j++)
        {
            stream << str.setNum(HMN[j]) << endl;
        }
        stream << "Checking: \n";
        gsl_vector_complex_set_zero(Check);
        for (int j=0; j<N; j++)
        {
            GSL_SET_COMPLEX(&z,0,0);
            for (int q=0; q<N; q++)
            {
                z=gsl_complex_add(z,gsl_complex_mul(gsl_matrix_complex_get(Acopy,j,q),gsl_vector_complex_get(X,q)));
            }
            gsl_vector_complex_set(Check,j,z);
            stream << str.setNum(GSL_REAL(gsl_vector_complex_get(Check,j))) << "\t" << str1.setNum(GSL_IMAG(gsl_vector_complex_get(Check,j))) << endl;
        }
    }
    writefile.close();
    gsl_permutation_free (P);
}

//! \fn Checking ortogonalisation
void Vainshtein::ortodox(QTextBrowser *stack)
{
    int N=100,m=0,n,coef;
    double r,ort[4],value;
    QString str,str1;
    for (int l=0; l<4; l++) ort[l]=0;
    for (int j=0; j<2; j++)
    {
        n=m+j;
        stack->append("m="+str1.setNum(m+1)+", n="+str.setNum(n+1));
        for (int i=0; i<=N; i++)
        {
            if(i==0 || i==N) coef=1;
            else if(i%2!=0) coef=4;
            else coef=2;
            r=a/N*i;
            ort[0]+=r*gsl_sf_bessel_J1(VMN[m]/a*r)*gsl_sf_bessel_J1(VMN[n]/a*r)*coef;
            ort[1]+=r*gsl_sf_bessel_J1(MMN[m]/a*r)*gsl_sf_bessel_J1(MMN[n]/a*r)*coef;
            r=a+(Shield-a)/N*i;
            value=gsl_sf_bessel_Jn(1,EMN[m]/a*r)*gsl_sf_bessel_Yn(1,EMN[m])-gsl_sf_bessel_Jn(1,EMN[m])*gsl_sf_bessel_Yn(1,EMN[m]/a*r);
            value*=gsl_sf_bessel_Jn(1,EMN[n]/a*r)*gsl_sf_bessel_Yn(1,EMN[n])-gsl_sf_bessel_Jn(1,EMN[n])*gsl_sf_bessel_Yn(1,EMN[n]/a*r);
            ort[2]+=r*value*coef;
            value=gsl_sf_bessel_Jn(1,HMN[m]*r)*gsl_sf_bessel_Yn(1,HMN[m]*a)-gsl_sf_bessel_Jn(1,HMN[m]*a)*gsl_sf_bessel_Yn(1,HMN[m]*r);
            value*=gsl_sf_bessel_Jn(1,HMN[n]*r)*gsl_sf_bessel_Yn(1,HMN[n]*a)-gsl_sf_bessel_Jn(1,HMN[n]*a)*gsl_sf_bessel_Yn(1,HMN[n]*r);
            ort[3]+=r*value*coef;
        }
        ort[0]*=a/N/3;    // E01
        ort[1]*=a/N/3;    // H01
        ort[2]*=(Shield-a)/N/3;   // coax E01
        ort[3]*=(Shield-a)/N/3;   // coax H01
        for (int l=0; l<4; l++) stack->append("I["+str1.setNum(l)+"]=\t"+str.setNum(ort[l]));
        if(j==0)
        {
            ort[0]=pow(a,2)/2*((1-pow(VMN[m],-2))*pow(gsl_sf_bessel_J1(VMN[m]),2)+pow(besselS(1,VMN[m]),2));
            ort[1]=pow(a,2)/2*pow(besselS(1,MMN[m]),2);
            double value1;
            value=gsl_sf_bessel_Jn(1,EMN[m]/a*Shield)*gsl_sf_bessel_Yn(1,EMN[m])-gsl_sf_bessel_Jn(1,EMN[m])*gsl_sf_bessel_Yn(1,EMN[m]/a*Shield);
            value1=besselS(1,EMN[m]/a*Shield)*gsl_sf_bessel_Yn(1,EMN[m])-gsl_sf_bessel_Jn(1,EMN[m])*besselSY(1,EMN[m]/a*Shield);
            ort[2]=pow(Shield,2)/2*((1-pow(EMN[m]/a*Shield,-2))*pow(value,2)+pow(value1,2));
            value=gsl_sf_bessel_Jn(1,EMN[m])*gsl_sf_bessel_Yn(1,EMN[m])-gsl_sf_bessel_Jn(1,EMN[m])*gsl_sf_bessel_Yn(1,EMN[m]);
            value1=besselS(1,EMN[m])*gsl_sf_bessel_Yn(1,EMN[m])-gsl_sf_bessel_Jn(1,EMN[m])*besselSY(1,EMN[m]);
            ort[2]-=pow(a,2)/2*((1-pow(EMN[m],-2))*pow(value,2)+pow(value1,2));
            value=besselS(1,HMN[m]*Shield)*gsl_sf_bessel_Yn(1,HMN[m]*a)-gsl_sf_bessel_Jn(1,HMN[m]*a)*besselSY(1,HMN[m]*Shield);
            ort[3]=pow(Shield,2)/2*pow(value,2);
            value=besselS(1,HMN[m]*a)*gsl_sf_bessel_Yn(1,HMN[m]*a)-gsl_sf_bessel_Jn(1,HMN[m]*a)*besselSY(1,HMN[m]*a);
            ort[3]-=pow(a,2)/2*pow(value,2);
            for (int l=0; l<4; l++) stack->append("I["+str1.setNum(l)+"]=\t"+str.setNum(ort[l]));
        }
    }
}

double Integrals::besselSY(int m, double x)
{
    return m/x*gsl_sf_bessel_Yn(m,x)-gsl_sf_bessel_Yn(m+1,x);
}

double Vainshtein::zerosCoaxE(int m, int n)
{
    int flag=0;
    double xi0=0;
    double xa,xb,det=0.05,val[2];
    val[1]=gsl_sf_bessel_Jn(m,(xi0+det)*Shield/a)*gsl_sf_bessel_Yn(m,(xi0+det))-gsl_sf_bessel_Jn(m,(xi0+det))*gsl_sf_bessel_Yn(m,(xi0+det)*Shield/a);
    do
    {
        xi0+=det;
        xa=xi0;
        val[0]=val[1];
        xb=xa+det;
        val[1]=gsl_sf_bessel_Jn(m,xb*Shield/a)*gsl_sf_bessel_Yn(m,xb)-gsl_sf_bessel_Jn(m,xb)*gsl_sf_bessel_Yn(m,xb*Shield/a);
        if (val[0]*val[1]<0)
        {
            flag++;
            if (flag==n)
            {
                xi0=(xa+xb)/2;
                for (int j=0;j<16;j++)
                {
                    val[0]=gsl_sf_bessel_Jn(m,xa*Shield/a)*gsl_sf_bessel_Yn(m,xa)-gsl_sf_bessel_Jn(m,xa)*gsl_sf_bessel_Yn(m,xa*Shield/a);
                    val[1]=gsl_sf_bessel_Jn(m,xi0*Shield/a)*gsl_sf_bessel_Yn(m,xi0)-gsl_sf_bessel_Jn(m,xi0)*gsl_sf_bessel_Yn(m,xi0*Shield/a);
                    if (val[0]*val[1]<0) xb=xi0;
                    else xa=xi0;
                    xi0=(xa+xb)/2;
                }
            }
        }
    }
    while (flag<n);
    return xi0;
}

double Vainshtein::zerosCoaxH(int m, int n)
{
    int flag=0;
    double xi0=0;
    double xa,xb,det=0.05,val[2];
    val[1]=besselS(m,(xi0+det))*besselSY(m,(xi0+det)*Shield/a)-besselS(m,(xi0+det)*Shield/a)*besselSY(m,(xi0+det));
    do
    {
        xi0+=det;
        xa=xi0;
        val[0]=val[1];
        xb=xa+det;
        val[1]=besselS(m,xb)*besselSY(m,xb/a*Shield)-besselS(m,xb/a*Shield)*besselSY(m,xb);
        if (val[0]*val[1]<0)
        {
            flag++;
            if (flag==n)
            {
                xi0=(xa+xb)/2;
                for (int j=0;j<16;j++)
                {
                    val[0]=besselS(m,xa)*besselSY(m,xa/a*Shield)-besselS(m,xa/a*Shield)*besselSY(m,xa);
                    val[1]=besselS(m,xi0)*besselSY(m,xi0/a*Shield)-besselS(m,xi0/a*Shield)*besselSY(m,xi0);
                    if (val[0]*val[1]<0) xb=xi0;
                    else xa=xi0;
                    xi0=(xa+xb)/2;
                }
            }
        }
    }
    while (flag<n);
    return xi0;
}

double Vainshtein::amq(double gamma, double alpha)
{
    int N=1000,coef;
    double Int=0,r;
    double x2m;
    for (int j=0; j<=N; j++)
    {
        if(j==0 || j==N) coef=1;
        else if(j%2!=0) coef=4;
        else coef=2;
        r=a/N*j;
        x2m=r*gsl_sf_bessel_J1(gamma/a*r)*gsl_sf_bessel_J1(alpha/Shield*r);
        Int+=x2m*coef;
    }
    Int=Int*a/N/3;
    return Int;
}

double Vainshtein::bnq(double chi, double alpha)
{
    int N=1000,coef;
    double Int=0,r;
    double FF,x2m;
    for (int j=0; j<=N; j++)
    {
        if(j==0 || j==N) coef=1;
        else if(j%2!=0) coef=4;
        else coef=2;
        r=a+(Shield-a)/N*j;
        FF=gsl_sf_bessel_J1(chi*r)*gsl_sf_bessel_Y1(chi*a)-gsl_sf_bessel_J1(chi*a)*gsl_sf_bessel_Y1(chi*r);
        x2m=r*gsl_sf_bessel_J1(alpha/Shield*r)*FF;
        Int+=x2m*coef;
    }
    Int=Int*(Shield-a)/N/3;
    return Int;
}

void viewinfo (QTextStream *stream, int N, gsl_matrix_complex *A, gsl_vector_complex *B, gsl_vector_complex *X, gsl_permutation *P)
{
    QString str,str1;
    *stream << "matrix A:\n";
    for (int j=0; j<N; j++)
    {
        for (int q=0; q<N; q++)
        {
            *stream << str.setNum(GSL_REAL(gsl_matrix_complex_get(A,j,q))) << "+ i*" << str1.setNum(GSL_IMAG(gsl_matrix_complex_get(A,j,q))) << "\t";
        }
        *stream << endl;
    }
    *stream << "vector B:\n";
    for (int j=0; j<N; j++)
    {
        *stream << str.setNum(GSL_REAL(gsl_vector_complex_get(B,j))) << "+ i*" << str1.setNum(GSL_IMAG(gsl_vector_complex_get(B,j))) << endl;
    }
    *stream << "matrix P:\n";
    for (int j=0; j<N; j++)
    {
        *stream << str.setNum(gsl_permutation_get(P,j)) << endl;
    }
    *stream << "vector X:\n";
    for (int j=0; j<N; j++)
    {
        *stream << str.setNum(GSL_REAL(gsl_vector_complex_get(X,j))) << "\t" << str1.setNum(GSL_IMAG(gsl_vector_complex_get(X,j))) << endl;
    }
}

//void Vainshtein::wenger2(QTextBrowser *stack)
//{
//    int N;
//    QString str,str1;
//    //QFile writefile("SOLVE.CSV");
//    //QTextStream stream(&writefile);
//    N=difQ[1]+difQ[2];    //without TEM
//    gsl_complex z,Res,delta;
//    double beta;
//    stack->append("k="+str.setNum(k));
//    stack->append("N="+str.setNum(N));
//    stack->append("II:\t"+str.setNum(pow(k,2)-pow(MMN[0]/a,2)));
//    stack->append("III:\t"+str.setNum(pow(k,2)-pow(HMN[0]/a,2)));
//    beta=sqrt(pow(k,2)-pow(MMN[0]/a,2));
//    GSL_SET_COMPLEX(&z,0,-beta/M_PI*2*(a*log(a/(Shield-a))-Shield*log(Shield/(Shield-a))));
//    Res=gsl_complex_exp(z);
//    Res=gsl_complex_mul_real(Res,-1);
//    for (int j=1; j<difQ[1]; j++)
//    {
//        GSL_SET_COMPLEX(&delta,sqrt(pow(MMN[j]/a,2)-pow(k,2)),beta);
//        Res=gsl_complex_mul(Res,gsl_complex_div_real(gsl_complex_pow_real(delta,2),pow(gsl_complex_abs(delta),2)));
//    }
//    for (int j=0; j<difQ[1]; j++)
//    {
//        GSL_SET_COMPLEX(&delta,0,-sqrt(pow(k,2)-pow(MMN[j]/Shield,2))-beta);
//        Res=gsl_complex_mul(Res,gsl_complex_div_real(gsl_complex_pow_real(delta,2),pow(gsl_complex_abs(delta),2)));
//    }
//    for (int j=0; j<difQ[2]; j++)
//    {
//        GSL_SET_COMPLEX(&delta,0,sqrt(pow(k,2)-pow(HMN[j],2))-beta);
//        Res=gsl_complex_mul(Res,gsl_complex_div_real(gsl_complex_pow_real(delta,2),pow(gsl_complex_abs(delta),2)));
//    }
//    stack->append("Coef=\t"+str1.setNum(GSL_REAL(Res))+"+i"+str.setNum(GSL_IMAG(Res)));
//    stack->append("|Coef|:\t"+str.setNum(gsl_complex_abs(Res)));
//    stack->append("arg(Coef):\t"+str.setNum(gsl_complex_arg(Res)));
//}

//! new method debug
//void Vainshtein::analitic(QTextBrowser *stack)
//{
//    int N;
//    QString str,str1;
//    QFile writefile("SOLVE.CSV");
//    QTextStream stream(&writefile);
//    //initialisation
//    N=difQ[1]+difQ[2];    //without TEM
//    gsl_complex z,Aq,beta1,beta,delta;
//    gsl_matrix_complex * A = gsl_matrix_complex_alloc (N, N);
//    gsl_vector_complex * B = gsl_vector_complex_alloc (N);
//    gsl_vector_complex * X = gsl_vector_complex_alloc (N);
//    gsl_vector_complex * Xnew = gsl_vector_complex_alloc (N);
//    int s;
//    gsl_permutation * P = gsl_permutation_alloc (N);
//    if(writefile.open(QIODevice::WriteOnly))
//    {
//        stack->append("Checking...");
//        for (int j=0; j<N; j++)
//        {
//            GSL_SET_COMPLEX(&delta,(double)j/44,1);
//            for (int q=0; q<difQ[1]; q++)
//            {
//                GSL_SET_COMPLEX(&beta,1,q+1);
//                gsl_matrix_complex_set(A,j,q,gsl_complex_pow_real(gsl_complex_sub(beta,delta),-1));
//            }
//            for (int q=0; q<difQ[2]; q++)
//            {
//                GSL_SET_COMPLEX(&beta,2,-q-1);
//                gsl_matrix_complex_set(A,j,difQ[1]+q,gsl_complex_pow_real(gsl_complex_sub(beta,delta),-1));
//            }
//            GSL_SET_COMPLEX(&beta1,1,1);
//            gsl_vector_complex_set(B,j,gsl_complex_pow_real(gsl_complex_add(beta1,delta),-1));
//        }
//        gsl_matrix_complex * Acopy = gsl_matrix_complex_alloc (N, N);
//        gsl_matrix_complex_memcpy(Acopy,A);
//        gsl_linalg_complex_LU_decomp (A, P, &s);
//        gsl_linalg_complex_LU_solve (A, P, B, X);
//        viewinfo (&stream, N, Acopy, B, X, P);
//        //! new method
//        for(int p=0; p<difQ[1]; p++)
//        {
//            GSL_SET_COMPLEX(&Aq,1,0);
//            GSL_SET_COMPLEX(&beta,1,p+1);
//            for(int j=0; j<N; j++)
//            {
//                GSL_SET_COMPLEX(&delta,(double)j/44,1);
//                z=gsl_complex_sub(beta,delta);
//                z=gsl_complex_div(z,gsl_complex_add(beta1,delta));
//                Aq=gsl_complex_mul(Aq,z);
//            }
//            for(int j=0; j<difQ[1]; j++)
//            {
//                if(j!=p) {
//                    GSL_SET_COMPLEX(&delta,1,j+1);
//                    z=gsl_complex_add(beta1,delta);
//                    z=gsl_complex_div(z,gsl_complex_sub(beta,delta));
//                    Aq=gsl_complex_mul(Aq,z); }
//            }
//            for(int j=0; j<difQ[2]; j++)
//            {
//                GSL_SET_COMPLEX(&delta,2,-j-1);
//                z=gsl_complex_add(beta1,delta);
//                z=gsl_complex_div(z,gsl_complex_sub(beta,delta));
//                Aq=gsl_complex_mul(Aq,z);
//            }
//            z=gsl_vector_complex_get(X,p);
//            stack->append(str.setNum(GSL_REAL(z))+"+i"+str1.setNum(GSL_IMAG(z)));
//            stack->append(str.setNum(GSL_REAL(Aq))+"+i"+str1.setNum(GSL_IMAG(Aq)));
//            gsl_vector_complex_set(Xnew,p,Aq);
//        }
//        for(int p=0; p<difQ[2]; p++)
//        {
//            GSL_SET_COMPLEX(&Aq,1,0);
//            GSL_SET_COMPLEX(&beta,2,-p-1);
//            for(int j=0; j<N; j++)
//            {
//                GSL_SET_COMPLEX(&delta,(double)j/44,1);
//                z=gsl_complex_sub(beta,delta);
//                z=gsl_complex_div(z,gsl_complex_add(beta1,delta));
//                Aq=gsl_complex_mul(Aq,z);
//            }
//            for(int j=0; j<difQ[1]; j++)
//            {
//                GSL_SET_COMPLEX(&delta,1,j+1);
//                z=gsl_complex_add(beta1,delta);
//                z=gsl_complex_div(z,gsl_complex_sub(beta,delta));
//                Aq=gsl_complex_mul(Aq,z);
//            }
//            for(int j=0; j<difQ[2]; j++)
//            {
//                if(j!=p) {
//                    GSL_SET_COMPLEX(&delta,2,-j-1);
//                    z=gsl_complex_add(beta1,delta);
//                    z=gsl_complex_div(z,gsl_complex_sub(beta,delta));
//                    Aq=gsl_complex_mul(Aq,z); }
//            }
//            z=gsl_vector_complex_get(X,p+difQ[1]);
//            stack->append(str.setNum(GSL_REAL(z))+"+i"+str1.setNum(GSL_IMAG(z)));
//            stack->append(str.setNum(GSL_REAL(Aq))+"+i"+str1.setNum(GSL_IMAG(Aq)));
//            gsl_vector_complex_set(Xnew,difQ[1]+p,Aq);
//        }
//        //! checking
//        stream << "Checking first: \n";
//        for (int j=0; j<N; j++)
//        {
//            GSL_SET_COMPLEX(&z,0,0);
//            for (int q=0; q<N; q++)
//            {
//                z=gsl_complex_add(z,gsl_complex_mul(gsl_matrix_complex_get(Acopy,j,q),gsl_vector_complex_get(X,q)));
//            }
//            stream << str.setNum(GSL_REAL(z)) << "\t" << str1.setNum(GSL_IMAG(z)) << endl;
//        }
//        stream << "Checking new: \n";
//        for (int j=0; j<N; j++)
//        {
//            GSL_SET_COMPLEX(&z,0,0);
//            for (int q=0; q<N; q++)
//            {
//                z=gsl_complex_add(z,gsl_complex_mul(gsl_matrix_complex_get(Acopy,j,q),gsl_vector_complex_get(Xnew,q)));
//            }
//            stream << str.setNum(GSL_REAL(z)) << "\t" << str1.setNum(GSL_IMAG(z)) << endl;
//        }
//    }
//    writefile.close();
//    gsl_permutation_free (P);
//}

